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Primordial magnetic fields and massive neutrinos can leave an interesting signal in the CMB 
temperature and polarization. We perform a systematic analysis of general perturbations in the 
radiation-dominated universe, accounting for any primordial magnetic field and including leading- 
order effects of the neutrino mass. We show that massive neutrinos qualitatively change the large- 
scale perturbations sourced by magnetic fields, but that the effect is much smaller than previously 
claimed. We calculate the CMB power spectra sourced by inhomogeneous primordial magnetic 
fields, from before and after neutrino decoupling, including scalar, vector and tensor modes, and 
consistently modelling the correlation between the density and anisotropic stress sources. In an 
appendix we present general series solutions for the possible regular primordial perturbations. 



I. INTRODUCTION 

The origin of the 10 _6 G magnetic fields observed in 
galaxies and clusters poses something of a problem for 
contemporary astrophysics pQ. Recent observations of 
galaxies at redshift z ~ 0.7 — 2 seem to show the fields 
were of comparable strength when the Universe was much 
younger, disfavouring a large dynamo amplification from 
tiny ~ 10 _20 G seed fields [3J [5] . Tentative observations 
of magnetic fields in elliptical galaxies, and a detection 
in a dwarf galaxy, also disfavour several dynamo mecha- 
nisms because they have little coherent rotation (see [4] 
and references within). There is also some evidence for 
~ 10~ 8 G fields coherent on megaparsec scales |S]. It 
may be possible to explain these observations in terms 
of astrophysically generated seed fields. Another inter- 
esting possibility is a primordial seed field. A primordial 
B ~ 10~ 9 G (comoving) field could lead to the observed 
galactic fields via adiabatic contraction alone, and might 
leave an interesting observable signature in the CMB. In 
this paper we revisit the calculation of the CMB power 
spectrum from primordial inhomogeneous magnetic fields 
and work towards robust theoretical predictions that can 
be used to test the primordial field scenario with CMB 
data. Since primordial magnetic fields are expected to be 
exponentially small in most early-universe models, any 
detection would be a clear signature of something very 
interesting. 

If a primordial inhomogeneous magnetic field is present 
it sources scalar, vector and tensor modes, giving rise to a 
signal in the CMB temperature as well as E and B-mode 
polarization. Previous calculations have indicated that 
~ 10~ 9 G fields (comoving) are detectable [51 [7], but have 
been incomplete in several respects. One complication is 
that the anisotropic stress due to the magnetic fields be- 
comes compensated by the neutrino anisotropic stress [5] , 
significantly reducing the perturbations sourced on large 
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scales after neutrino decoupling. Recent work by Kojima 
et al. [9] has claimed that the presence of massive neu- 
trinos leads to a significant change in this compensation 
mechanism, giving rise to a dramatic enhancement of up 
to eight orders of magnitude on the large-scale E-mode 
polarization power spectrum. For interesting neutrino 
masses this would, if true, be a clear signal of primordial 
magnetic fields. Clearly this claim merits further inves- 
tigation, though we shall ultimately show that the ef- 
fect, though interesting, is much smaller than previously 
claimed. 

In the early universe the massive neutrinos are ex- 
pected to be relativistic, with the most massive eigen- 
state only becoming non-relativistic around recombina- 
tion or later [ID]. We perform a systematic analysis of 
the primordial perturbations to lowest order in the mass, 
generalizing previous results for the general primordial 
perturbation to the realistic case where one or more of 
the neutrinos is massive. This also allows us to calculate 
the series solutions consistently in the presence of mag- 
netic fields, and see the leading corrections due to the 
neutrino mass effect. We also discuss the tight coupling 
approximation, which is useful after the modes come in- 
side the horizon but before Thomson scattering becomes 
ineffective. 

The calculation of the CMB power spectrum from pri- 
mordial magnetic fields is further complicated because 
the scalar, vector and tensor sources are all quadratic 
in the underlying magnetic field, and the scalar modes 
have more than one source term. These sources are cor- 
related, so we show how to calculate the various source 
power spectra and correlations from the power spectrum 
of the magnetic field, and how to use these for a numer- 
ical calculation. 

In this work we address the effects of magnetic fields 
in sourcing primary anisotropies in the CMB, however, 
magnetic fields present after recombination also have an 
observational effect on the polarisation by inducing Fara- 
day rotation [TTJ [T^]. Such rotation converts E-mode 
polarisation into B-modes with a strong dependence on 
the frequency of the radiation oc Bjv 2 but it is small at 
the usual frequencies for CMB observation, and we will 
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neglect it in our analysis. 

Throughout this work we use a 3+1 splitting of General 
Relativity, working with a gauge invariant linear pertur- 
bation theory similar to that of Bardeen |13j and Durrer 
[IH El] • Our choice of gauge invariant variables is chosen 
as a close analogy to the Conformal Newtonian Gauge 
(CNG). We use a metric 

ds 2 = a(r) 2 -(l + 2A)dr 2 -2B z drdx l 

+ (S i:j + 2Hij) dx i dx j (1) 

where we can further decompose Bi and into their 
scalar, vector and tensor contributions. Our decompo- 
sitions are performed in the same manner as Ref. |16j . 
In fc-space the decomposition for the rank-1 and rank-2 
three-tensors are written as (using Bi and Hij as exam- 
ples) 

Bi = BQf ] + B™Q™ 

= H L 8ijQ^ + HxQij + H m Qf? + (2 >Q< 2) (2) 

where the harmonic Q {m) functions give the form of each 
perturbation type for a specific k mode, with m = giv- 
ing scalar perturbations, and m = 1, m = 2, vector and 
tensor respectively. In the above it should be understood 
that we implicitly sum over the two vector and two tensor 
modes, for example 



(3) 



whilst a quantity like -ff (1) appearing on its own can stand 
for either H i+1) or H { ~ 1) consistent with the context. The 
scalar harmonic functions are 



For reference the self contraction of this is 
Qij 2) *Q* 3 ' (±2) = §■ As we would expect, quantities 
of different types are always orthogonal, for example 

Q(±2)*Q ii(±1)=0 _ 

Generally we will drop the superscript on scalar per- 
turbations like X (a) , in favour of simply X. 

As a further illustration, let us examine the pertur- 
bations to the energy- momentum tensor T^ u . In the 
Conformal Newtonian Gauge, the gauge invariant quan- 
tities we use are exactly equivalent to perturbations of 
the Energy-Momentum tensor. The density A, velocity 
Vi, pressure it, and anisotropic stress IK perturbations 
are defined by 



T ° = -p(l + A) , 
T? = (p + p)V l7 



(8) 
(9) 
(10) 



where p and p are the density and pressure respectively. 
As above we will further decompose the three- vector and 
tensor quantities into the different perturbation types. 
The velocity three- vector decomposes as 



Vi = VQ?> + n m Q 



(11) 



where the vorticity is the vector-type velocity. The 
traceless anisotropic stress tensor becomes 



m = no^ + n (1 »Q^ + n (2) Q 



(12) 



Generally we will rewrite the pressure perturbation 7r in 
terms of an entropy type perturbation Y and the density 
A 



Q(O) = gZk-X 

Qfj = [fc- 2 v 4 v, + ^/3] g (0) 

The vector harmonics are 

Q 



(±i) _ (±) ik-x 



1, 



(±1) _ 



7'L-P (±) P ik ' X 



(4) 



(5) 



r + ^A. 

w 



(13) 



Although we will not explicitly demonstrate it, T is 
gauge- invariant. Above we have used w and c 2 , defined 
as w = p/p and the sound speed c 2 = p/p. 

We will restrict ourselves to a flat geometry throughout 
this work. 



II. NEUTRINO PERTURBATIONS 



where we decompose our vectors with the helicity basis 



A. Kinetic Theory 



V2 



(e, l ±<e?) 



(6) 



with e 1 and e 2 being unit vectors orthogonal to k. Note 
that e (±) -e (±) = e (±) * -e^) = 0, whilst e (± '* -e (±) = -e (±) • 



e (T) = i. F rom t his we can see that Q\f 1)m Q ij (±y = \. 
For the tensors we make the further definition of e 
\/3/2 e^ ±) e'. ±:) . Using this, the sole tensor harmonic is 



(±2) 



QT = ^ 2>e4k ' x 



To describe the behaviour of neutrinos in the early uni- 
verse, we must turn to the full machinery of the Boltz- 
mann equation. We start with the phase space distri- 
bution of the particle density on a spatial hypersurface, 
defined by 



dN = f v (x\Pj,T) d 3 xd 3 P 



(14) 



where Pi is the canonical 3-momentum, the spatial part 
(7) of the covariant 4-momentum P^. The primary quantity 
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we will require is the energy-momentum tensor which is 
determined from the distribution function f u by 



T"„ = 



d 3 P(-g) 

Pn 



1/2 



(15) 



Following the convention in the literature (HO [T7] we 
will rc-express the distribution function in terms of quan- 
tities in the frame of a comoving observer. We use the 
locally Minkowski tetrad satisfying g^ v ^a e b = Vab- In 
terms of the co-ordinate basis, and where we have avoided 
fixing a gauge 



e^a- 1 [(1-A)d -B i d i \ 

-l 



a 



(l-H L )d i -Hid j 



(16a) 



with Hf containing the trace-free scalar, vector and ten- 
sor contributions. This allows us to write the momentum 
in terms of quantities measured in the comoving tetrad 
P = P^d^ = Tr a e a , where ir° is the observed energy and 
7T 1 the momentum in that frame. Applying Hamilton's 
equations to the system implies that the conjugate mo- 
menta will remain constant in a purely FRW universe. 
This means the proper 4-momenta will decay away with 
a -1 . In order to remove this redshifting of the energy 
and momenta we will write them in terms of the scaled 
quantities e and q defined by 



7T 
7T' 



e/a , 
qri 1 1 a 



(17) 



where n l is the unit vector in the direction of the mo- 
mentum. Both e and q are constant on the background 
by definition. By considering P ■ P = Tr a n a we find a 
slightly modified energy momentum relation 



e(q) = (q 2 + a 2 m 2 )^ 2 



(18) 



Prior to their decoupling, neutrinos are in approximate 
thermal equilibrium with the rest of the Universe. Con- 
sidering only the unperturbed case for the moment, the 
phase-space distribution function of the neutrinos /„o will 
be Fermi-Dirac at a universal temperature. We expect 
this distribution to be isotropic and homogenous, and 
thus only be a function of the momentum magnitude q 
(in the guise of the comoving energy) and the time r (by 
virtue of the temperature). Therefore it takes the form 



Uo( a , T ) 



9, 



1 



^ e E(q)/k B T(7 



(19) 



where the neutrino energy measured by a comoving ob- 
server is E = e/a. As the temperature decreases with 
1/a, the combination E(q)/kBT = e(q)/kBTo is constant, 
depending on To, the temperature today. 

At neutrino decoupling, this distribution becomes 
frozen in. The neutrino mass is insignificant compared to 
any thermal energy, so its contribution can be neglected 



in the distribution function. This allows us to set e = q 
(within the distribution only), leaving the unperturbed 
function as 



1 



h^ p e q/k B T + 1 



(20) 



We will define the first order perturbations to the dis- 
tribution ijj v by 



fu(x\Pj,T) 



")] 



(21) 



with ip v containing the scalar, vector and tensor contri- 
butions. This quantity is gauge dependent; later we will 
form a gauge invariant equivalent. 



We want to rewrite the integral (15 1 in terms of our 



comoving quantities, retaining terms up to first order. 
Firstly, the term d 3 P (—gY^/Po forms a co-ordinate in- 
variant measure for the integration, and can be re- written 
in terms of the comoving quantities 



-nW 2 



d 3 P(-g) 
~P7 



— a 2 dqdQ Tr 



(22) 



This removes the metric perturbations contained within 
the integration measure. Re-expressing the P M P„ gen- 
erates a plethora of terms, including terms first-order in 
the metric perturbations. However these terms all de- 
pend upon a single power of the momentum direction 
and couple only with the isotropic distribution /„o; they 
are thus eliminated by their symmetry. The remaining 
terms are simply 



P»P V = a - 2 (eSi; + qn l 5?)(-e5l + qn,^ v ) + 



(23) 



Decomposing into distinct components this leaves us with 



q 2 dqd£l n e f„o{q) [1 + ip u ] 



T" = a 4 / q 2 dqdtt n qrii f v ^ v , 



(24) 



TV = a' 4 I q 2 dqdn n ^-n l n 3 f v0 (q) [1 + 



valid in all gauges up to first order. 

The evolution of the distribution function is governed 
by the collisionless Boltzmann (or Vlasov) equation, 
which simply expresses that, without collisions, the num- 
ber of particles is conserved along a trajectory in phase 
space 



Df 
Dt 



df 
8t 



dx l df 
dr dx l 



dqdl 
dr dq 



drij df 
dr drij 



(25) 



In order to address the evolution of the perturbations to 
the distribution function, we need to separate this out 
into equations for the background, and the separate per- 
turbation types. We address this in subsequent sections. 
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B. Background Quantities 

At zeroth order the collisionless Boltzmann equation 
simply shows that the distribution remains constant, that 
is, f u o is independent of time. At this order the energy- 
momentum tensor can be described fully in terms of den- 
sity and the pressure. The density is given by 



p v = 4ira 



q 2 dqef v0 (q) 



(26) 



With a non-zero neutrino mass the pressure is no longer 
simply related to the density. It's instead 



Instead of the scalar metric perturbations we will use the 
gauge-invariant Bardeen potentials 



# = A 



B + HB) - 



) k 2 ( 



Ht + TLHt 



= A-Uajk- &/k 

* = -H L - ^H T - jB 
= -K + Ha/k . 



H 

k 2 



(31) 
(32) 



The potential \& should not be confused with the distri- 
bution perturbation 5'< / m) . Usually the context will make 
this clear. In the above, 7Z is the 3-Ricci scalar 



4?r 



Pv 



q dq—Uo(q) 



(27) 



As we would expect, the equation of state w v = p v / ' p v is 
still defined by the ratio of these two quantities, yielding 
a mass dependent w ^ 1/3. 



U = H L + ^H T 



(33) 



Written in terms of gauge-invariant quantities, the Boltz- 
mann equation becomes 



ikp-^ 



d In f v0 
dhiq 



$ - ikp-^ 

q 







(34) 



C. Scalar Perturbations 

Initially we will just address scalar perturbations to 
the distribution function i/^ 0) , considering vectors and 
tensors later on. 

We will perform a harmonic expansion of all our quan- 
tities. In flat-space this is just the Fourier transform. The 
Boltzmann equation ( 25 ) can be expanded at first order 



including all the (non-gauge fixed) metric perturbations 
[Hj giving 



dhiq 



iku-A + Bkfj 2 

q 



H L - (fi 2 - 1/3) H T 



(28) 



where /i = n % ki and the dot is the derivative with respect 
to conformal time r. To move to a gauge-invariant for- 
malism we follow Durrer and Straumann 14 and define 
a new gauge-invariant distribution perturbation 



■ 

ding 



H .e 
k q 



(29) 



where a is the shear on spatial hypersurfaces a — Ht /k— 
B, and % = a/a is the conformal Hubble parameter. 
This differs from Durrer's definition in that we have cho- 
sen our definition to coincide with the CNG result (added 
terms vanish in a zero-shear gauge). For comparison 
Durrer's invariant perturbation T and our definition are 
linked via 



(0) _ 



^ m _dlnf u0 



dlnq 



4> 



(30) 



As mentioned, our choice of gauge invariant variables is 
designed such that this is equivalent to the CNG version. 

The dependence on the momentum direction within 
the Boltzmann equation makes a direct solution tricky. 
We take the standard approach and expand out into an 
angular basis. Whilst for scalar perturbations it suffices 
to expand in the Legendre polynomials Pi (fi) , for vector 
and tensor perturbations it is much more convenient to 
use a method similar to Ref. |16j . where we expand out 
into spherical harmonics YJ m (#,0). Under this expan- 
sion, the different types of perturbations are separated 
by their m value, with scalar (m = 0), vector (m = 1) 
and tensor (m = 2) modes all evolving separately in the 
usual manner. Expanding out the entire distribution per- 
turbation 

oo / I ~. 

*" = £ EH)V2iTi^ )(fci '^ m(ni) ' (35) 

1=0 m=0 

we can relate the momentum integrals of the multipole 
moments to the standard gauge invariant perturbations 



4?r 
p v a 



, , q 2 dqeU(q)K°o(h,q) , 



47T 



, , q 2 dqqf uQ (q)^l(k h q) , (36) 



4-7T 

5p u a 4 



q 2 dq^f u0 (q)^(k h q) 



D. Thermal Perturbations 

The most natural perturbation that could be set up is 
from a purely thermal distribution, where we perturb the 
neutrinos by having a position and direction dependent 
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change to the temperature. To take this into account 
let us re- write our distribution perturbation in a slightly 
different manner. The total distribution 

9s ' (37) 



f v (x\q,nj 



/j3 e q/k B T (i+e v 



1 



fuo(q) 



dlnfvo . i 

1 — 9„{x ,q,nj,T) 

dmq 



For generality we leave 9 a function of q. For a pure 
temperature perturbation it must be temperature inde- 
pendent. Relating this to our previous gauge invariant 
perturbation ty v we find that 



d In f v0 
ding 



9 u (x\q,n j ,T) + a 



H .e 

T +t-p 
k q 



(38) 

and from this construct a temperature-like gauge- 
invariant perturbation 

"H 



9„(a;\ 9,^,7-) = Q v + , 



— + i-p 
k q 



(39) 



such that = - d dllq Bv Substitution into the Boltz- 



mann equation (|34| produces 



ikp,-^ 



o 



(40) 



When 77i = or r— > 0, e — q and the Boltzmann 
equation becomes momentum independent; the pertur- 
bation remains purely thermal. However even perturba- 
tions which start in a purely thermal state evolve away 
from it as the mass becomes important. For this reason 
we must keep 8„ a function of q, though we will restrict 
ourselves to purely thermal initial conditions. 



III. MASS EXPANSION 

To treat massive neutrinos in the early Universe, when 
the mass is small in comparison to the typical momen- 
tum (approximately fc^X 1 ,,), we will expand the system to 
first order in the neutrino mass squared. This will allow 
us to directly tackle the integrated distribution function, 
making it possible to find initial conditions up to this or- 
der in the neutrino mass. For a more general approach 
see P) . 

In both the integrals for the energy-momentum tensor 
and the Boltzmann equation itself the mass dependence 
comes in from factors of e/q or its inverse. Expanding 
this out gives e/q = 1 + m 2 /2q 2 + ■ ■ ■ (with a minus for 
the inverse). For the background quantities performing 
this expansion gives 



p u = Aira 4 / q 2 dqej v0 (q) 
= 4ira~ i j ' q 2 dqqf uQ (q) 

{-, 1 -2 2 \ 
= PvO 1 + ~ m a + ' 



1 ma 



2^2 



2 q 2 



(41) 



where p v § is the density for massless neutrinos, and the 
scaled mass fh 2 = m 2 /q 2 , with the q 2 factor being defined 
via 



J_ = Jq 2 dqqq 2 /„o(g) 
q 2 fq 2 dqqf» a (q) 



(42) 



which is essentially the momentum averaged inverse 
square momentum. This depends only on the background 
distribution, and is time independent. Thus we factor it 
into the new dimcnsionless mass fh. In terms of back- 
ground quantities 



m 



10 m 
7tt 2 k B T 



(43) 



As with the density we expand up to m 2 for the back- 
ground pressure p„ and the equation of state w v giving 
the form 



i'r = I'm ( 1 - 7, to2 ° 2 



w v = Pu/Pu = - (1 - m 2 a 2 ) 



(44) 



The perturbed quantities are slightly more tricky. For 
example taking A„ and expanding out in mass gives, 



L±V — 

p v a 



1 TV^C? ^ 

i , <?dqqf M {q)y v0 [l + --^) , (45) 



with a similar pattern for the other perturbations. 
Schematically we have two types of terms we need to 
integrate; J q 2 dqqf v ^ v and J qdqf U Q^ v . For a moment 
let us consider what happens to these integrals for ther- 
mal perturbations in the case of massless neutrinos 



q 2 dqqf v0 ^ u = -Q v I q 2 dqqf, 



dhiq 



49, 



a 4 PuO 
An 



and similarly the second integral equates to 



qdqf v0 ^ v = 29„ 



4 \ 1 

a p v g \ 1 

47T 



(46) 



(47) 



From this we can make the connection that, at zeroth 
order in the mass expansion, the second integral is linked 
to the first via 



1 

2<f 



qdqfvoVu = 7^2 I <fdqqf v0 ^ u + 0(fh 2 ) , (48) 



and as the second integral appears at first order in m 2 in 
the mass expansion, we can use this relation to simplify 
the expression for A„ above, giving 



A„ = -^j ( I + >/V ) / l fd qq j rt] *„„ . (4!)) 
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This happens similarly with the other perturbed quanti- 
ties, and allows us to follow the convention of forming a 
momentum integrated function 



F(ki,/j,,T) = 



Jq 2 dqqf v o(q)^ v (ki,q,n) 



J q 2 dqqf„ (q) 



q 2 dqqfuo{q)^u(ki,q,fJ,) . (50) 



As with the distribution perturbation, we will expand F 
into spherical harmonics 



F(k U f,,T)=f2(-i) l J^jFr(k i ,T)Y^ t i). (51) 

1=0 

Each moment F\ takes the same form as the integrated 
F of (50 1 with the being replaced by ^ v i. 



These lead to a succinct form for the perturbations 



For the dipole (I = 1) 

Fi + | (l - ^m 2 aA [2F 2 - 5F ] = (4 + m 2 a 2 )k^ . 

(56b) 

Finally for the quadrupole and higher moments (I > 2) 



Fi+kll 



' ^.2 2 

-m a 



l + l 



21 



F, 



I 



i+i 



21 - 1 



Fi-i 



= . 

(56c) 

As should be expected sending fh — > takes everything 
to the well known massless case. 

Taking this mass expansion to higher order becomes 
more difficult: Taylor expanding the background quanti- 
ties inside the integral produces divergent integrals at or- 
der m 4 and above. We show how to do the higher-order 
expansion in Appendix |Aj this shows that the leading- 
order mass expansions are correct to 0(m 4 log(m)). 



PuqFq (l + \m 2 a 2 
Po (l + \m 2 a 2 ) 



= F„ ( I- -in 2 a 2 



(52) 



to order m 2 . Similarly, 



1 



n„ = -f 2 

5 



i x - 2 2 

1 ma 

4 



1 



- 2 2 

-m a 



(53) 



All that we need to do to have a working prescription for 
calculating the massive neutrino evolution is to turn the 
Boltzmann equation into a hierarchy for solving for the 
Fi . First we take the Boltzmann equation and expand 
to first order in m 2 , then integrate it over J q 2 dq q f v o(q) 
and divide by the same quantity to produce an equation 
for the evolution of F. We employ the same trick as above 
(in Eq. ( [48] )) to turn the m 2 /q 2 quantities into terms in 
F. This results in 



F 



ikfiF 



1 



- 2 2 

m a 



4<f> - 4ifc/j* [ 1 + ~^rn 2 a 2 



(54) 



We then substitute a spherical harmonic expansion for 
F, and using the identity that 




4(/ + 1) 2 -1 



(55) 



we obtain the hierarchies for F. Separating these out into 
coupled equations for each I gives three distinct cases. 
For the monopole (I — 0) 



IV. VECTOR PERTURBATIONS 

To consider vector perturbations we proceed down a 
similar line to the scalar perturbations. Though not en- 
tirely free of gauge issues, many of the complexities will 
disappear. Firstly, whilst there are two vector-type met- 
ric perturbations B {1) and H {1 \ we have one degree of 
gauge freedom, and so only one perturbation can be rel- 
evant. As before, rather than fixing a gauge we form 
one gauge-invariant variable. For vector perturbations, 
the shear-like perturbation er (1) = H lr> /k — B {1) is gauge- 
invariant and we use this as our metric variable. 

The vector contribution to the distribution function 
4'< / 1) is itself gauge-invariant, (see [H]), and thus the 
Boltzmann equation governing it is 



*£> +*fcji-*£ ) - d ^ 1 n i n^k i efka w 
e dmq J 







(57) 



There are two contributions to the energy-momentum 
tensor: the velocity, , which is gauge dependent, and 
the anisotropic stress n^ 1 ' which is gauge-invariant. As 
a gauge-invariant velocity, we use the neutrino vorticity 
Slj, 1 ' = w^ 1 ' — B {1) which is conveniently related to the dis- 
tribution perturbation. The two neutrino perturbations 
are therefore 



4tt 



3(p u +p„)a 4 



q 2 dqqf u v(q)^i(k l ,q) 



15p„a 4 



q 2 dq^-f v0 (q)^(ki,q) . (58) 



Performing the same momentum integration as for the 
scalars (with the same restriction to thermal modes), 
we find an equation governing the momentum-integrated 



F + -fi [ I 



1 



m 2 a 2 ) = 4$ 



(56a) 



F' 



ikfiF 



(i) 



i 1-22 

1 771 a 



-Y^ka w 



(59) 
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where we have used that in l n 3 kie (1) = y/An/lbY^. In- 
serting the spherical harmonic expansion, and using the 
identity (55 1, the moments of the Boltzmann equation 



are for I = 1 



*f > + (l ~ |m a o a ) = 



(60) 



for I = 2 



F^+kll- -m 2 a 2 



V3 



and for I > 2 



(61) 



A; ( 1 to 2 a 2 

4 



2/ 



2 _ 1 ,/Z 2 _ I 



2^ - 1 



z-1 



= . (62) 



are no relevant contributions from the I = 0, I = 1 and 
to 7^ ±2 terms in the sum. For I = 2 we have 



and for I > 2 we require 
1 



1 



- 2 2 

to a 



4# ( 



k ( 1 m 2 a 2 

4 



2Z + 3 
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The only tensor contribution to the energy momentum 
tensor comes from the anisotropic stress II„. This is eas- 
ily expressed in terms of the expanded F {2) with 



n£> = !(i + imV) r.r ! 



(69) 



Finally we need to rewrite both n^' and in terms of 
the integrated functions. These are nearly identical 
to the scalar equivalents, only with different coefficients 



nw(fci) 
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(63) 



V. TENSOR MODES 

As is well known, tensor perturbations are manifestly 
gauge invariant and so we need not concern ourselves 
with any gauge issues. Other than that we follow the 
same track as for the scalar perturbations. The Boltz- 
mann equation for tensor modes takes on the form 



e ding 



nWHlf = . (64) 



We then momentum integrate the equation to produce a 
single equation in terms of the F <2) . Again we have re- 
stricted ourselves to initially thermal perturbations, giv- 



ing 



F (2) +i 



k^F {2) ( 1 - -to 2 a 2 ) = -4nVe<f ] H {2) . (65) 



Using the helicity basis, the quantity n*n J e^ 3) is simply 
written in terms of spherical harmonics as 



-Y4 



(66) 



To finish off, we need to rewrite the Boltzmann equation 



(65) with a spherical harmonic decomposition. There 



VI. EVOLUTION EQUATIONS 

The behaviour of the early Universe is accurately de- 
scribed by linear perturbation theory, reducing to a sys- 
tem of coupled linear differential equations. We have 
discussed the perturbation equations for the neutrinos in 
previous sections. Here we briefly describe the remaining 
equations for the evolution of the metric potentials and 
the other matter species. 

The 3 + 1 splitting of the Einstein equation = 
&TtGT^ v decomposes into sets of equations for each of the 
scalar, vector and tensor contributions. For the scalar 
perturbations we have four equations generated by the 
splittings. There are two equations formed by the (00) 
and (Oi) components 



fc 2 $ = 



A + 3(l + w)^y 
k 



k($ + Hy) = -H 2 {l + w)V 



(70a) 
(70b) 



where the first is the equivalent of the classical Poisson 
equation. The spatial part (ij) splits into two further 
equations from the trace and traceless parts. The equa- 
tion from the trace is 

$ + W(¥ + 2$) + (m + -h 2 ) + ifc 2 ($ - #) 

= ^ 2 (c 2 A + w r), (70c) 

where T is the perturbation to the entropy of the system, 
and c 2 = p/p is the total sound speed of all the matter 
species. The final equation is from the traceless part 



fc 2 ($ - = m 2 wn. 



(70d) 



There are two vector equations, one from the (Oi) part, 
and the second from the vector contribution to the (ij) 
components: 



k{& ( 



k z a m = -m 2 (1 
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(71a) 
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same manner as for the neutrinos, with the distinction 
that they are are massless bosons, and interact with the 
baryons via Thomson scattering (see [23 E3). ^ e 
calculation requires a consistent treatment of polariza- 
tion; we do not repeat this here, see e.g. Ref. [TS] for the 
details. The photon hierarchy is concisely written as 



There is a single equation for the tensor modes 

H {2) + 2UR {2) + k 2 H {2) = m 2 wll {2) . (72) 

The matter evolution equations are well known and 
are most generally derived from the Boltzmann equation, 
here we will just give the results. We consider the stan- 
dard three matter species beyond neutrinos: baryons, 
photons and cold dark matter, giving their perturbations 
in terms of A, V and II as before. 

The matter species have essentially no velocity disper- 
sion, the anisotropic stress and higher momentum mo- 
ments are all zero. Hence they contribute only to scalar 
and vector modes and can be described entirely in terms 
of A, V and il. Simplest is dark matter as it has no 
interactions. For the scalars 



A c = -kV c + 3$ , 
V c = -HV c + kV , 



(73a) 
(73b) 
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the source terms S[ m} describe the interactions with the 
gravitational potentials and other matter species. The 
non-zero terms arc for the scalars 
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for the vectors 

S[ r> 

and for the tensors 
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ka {1 > +t~ 1 P< 1 \ (76b) 
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and for the vectors 



(73c) 



We can see that any vector solution for CDM must be 
decaying, and so we will neglect it. 

The baryons couple to the photons via Thomson scat- 
tering, but also interact with any magnetic field via the 
Lorentz force giving an extra source term 



A 6 = -kV b + 3<S> , 
V b = -HV b + kc 2 sb A b 



(74a) 



■ k* + RT~ l (Vy - Vb) 



-kR 



(74b) 



where P (m) is the anisotropic Thomson source and con- 
tains the coupling to the polarisation 



l r 

10 



(77) 



In terms of the photon multipole moments the usual mat- 
ter sources are 



A 7 = 9 { q\ 



n (1) 
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(78) 



where the baryon sound speed is c 2 b = Sp b /Sp b . The two 
quantities A g and 11^ are the magnetic equivalents of the 
density and anisotropic stress perturbations. We make 
a thorough definition in the next section. The vector 
equation is 



■Rr^(Sl^-n^)--Ru>^, (74c) 



where R = 4 ( o 7 /3pb and r c is the timescale for Thom- 
son scattering, the inverse of the opacity, t^ 1 — an e aT- 
As with CDM there are no tensor perturbations to the 
baryon distribution. 

Describing the photon perturbations requires the full 
mechanics of the Boltzmann distribution. Constructing 
the gauge invariant perturbation equations is done in the 



A. Regular initial conditions 

We use the full system of perturbation equations to 
calculate initial series solutions for modes well outside 
the horizon in the early radiation-dominated epoch, af- 
ter neutrino decoupling but well before recombination. 
These are needed to provide the correct initial conditions 
for Boltzmann codes such as Camb [T!5] and Cmbfast 
|20) . We have calculated the complete set of all known 
regular modes for the standard matter species (dark mat- 
ter, baryons, photons and neutrinos) for the scalar, vector 
and tensor type perturbations. We have also included all 
the compensated magnetic modes for three perturbations 
types. By using our expanded neutrino equations (see 
Section III) we can include neutrinos of non- negligible 
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mass, with solutions accurate to order m 2 . Thus our so- 
lutions include both massless neutrinos and a number of 
degenerate massive species. 

We make several standard approximations, firstly we 
assume we are in the regime of tight coupling between 
photons and baryons where Thomson scattering prevents 
slippage between the fluids, giving V b w Vy (see [TTjV 
This gives two parameters which much be small: there 
must be many scatterings per wavelength of the pertur- 
bation kr c <C 1; and the scattering rate must be large 
compared to the expansion rate r c /r <C 1. We take the 
leading order corrections to this and truncate the tight 
coupling hierarchy by assuming the photon anisotropic 
stress IT 7 is negligible (it is suppressed by a factor fcr c rel- 
ative to the velocity). We also assume that the baryons 
are pressureless with Wb = c 2 b = 0, neglect any change 
in the background ionization fraction and degrees of free- 
dom, and as before assume a flat universe. Standard dark 
energy does not affect the result until C(r 5 ). 

The solutions are too lengthy to list in the main text 
and so we include them in Appendix [BJ 



VII. PRIMORDIAL MAGNETIC FIELDS 

We will consider a stochastic background of magnetic 
fields B l (x J , t) generated by some mechanism in the very 
early Universe. As for all the periods of interest the Uni- 
verse contains a highly ionized plasma, Maxwell's equa- 
tions at first order show that the field is frozen in, with an 
amplitude decaying with 1 /a 2 . From this we separate out 
the time evolution and write B 1 (x : ',t) = B l (x^) / o,(t) 2 . 
For a thorough discussion of the dynamics of cosmolog- 
ical magnetic fields, see [21]. The non-zero components 
of the energy-momentum tensor are 



rpV 

J- n — : 



1 



r^ 2 (x) , 



(79a) 



As there is no magnetic field on the background, the per- 
turbations of the stochastic background are manifestly 
gauge invariant. We construct two perturbations As and 
IIb, defined by 



1 o 



(80a) 



where we include the factors of p 7 and p 7 to take account 
of the a -4 factors. As usual II* can be decomposed in 
the standard manner into scalar, vector and tensor con- 
tributions. 



A. Magnetic Modes 

Though the exact mechanism by which magnetic fields 
may be produced in the primordial Universe is unclear, 



we are still able to address their observational conse- 
quences. We imagine that the production of magnetic 
fields occurs quickly at some time tb, prior to the de- 
coupling of neutrinos from the photons at time t„. We 
assume that this decoupling is effectively instantaneous. 
Below we briefly review what happens for the scalar case. 
This is discussed in detail in |22| using the synchronous 
gauge, where the calculations are somewhat simpler. Our 
gauge-invariant notation has the difficulty that some of 
the perturbations diverge on the superhorizon scales we 
are interested in, and this needs to be carefully addressed. 
The Mathematica notebook used for the calculations of 
the gauge-invariant scalar, and tensor case can be found 
at |http://camb.~in fo/jrs/, and describes these issues 
in more detail. 



Combining the four scalar Einstein equations ( 70 1 al- 



lows us to form the Bardeen equation for the potential 
$ which is sourced only by the total anisotropic stress n 
and the entropy T 

$ + 3K(1 + c 2 )$ + [3(c 2 - w)H 2 + c 2 fc 2 ]$ 



= 3w 



k 2 



k 2 



r 



k 2 

■nu- — n 

■ > 



-2-Hn + 3-H 2 (l-c 2 /w)n . (81) 



Prior to neutrino decoupling the Universe is dominated 
by the combined radiative fluid with c 2 — w = | . In this 
limit the Hubble parameter T~L = r _1 . The fluid is tightly 
bound to the trace amount of Baryons and cannot de- 
velop any anisotropic stress, and so the only anisotropic 
stress comes from the primordial magnetic source, the 
constant n b ■ Until neutrino decoupling there is no mech- 
anism to compensate this, and it will act as a source for 
the potentials. We will only discuss the anisotropic stress 
as the magnetic density perturbation must be compen- 
sated at generation on energy conservation grounds [23] . 
We reduce the Bardeen equation to the radiation domi- 
nated limit 



3fc 2 



r 2 $ 



4t$ 



+ fc 4 T 4 $ = -HyUn (6 + k 2 



(82) 



This can be solved exactly, and in the superhorizon limit 
of small kr it reduces to a solution of 



$(r) 



k 2 T 2 fc 3 T 3 



6kr 



c 2 --i? 7 n B log(r) (83) 



which has a singularity for kr — 0. As we are concerned 
with superhorizon modes, we check the physicality of this 
by examining the co-moving curvature perturbation £ = 
$ + 2(* + $/"H)/3(l + w) finding that 



C(r) = C(tb) 



7 llB 



log (t/tb) + 



TB 

2t 



(84) 



where we have absorbed the remaining constant terms 
by demanding continuity of the £ and the comoving den- 
sity perturbation (equivalent to continuity of $). All the 
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primordial contributions to the curvature are contained 
within C( t b)- 

At time r„ the neutrinos decouple from the radiative 
fluid. By considering their Boltzmann hierarchy we can 
examine what happens next. Combining the I = 1 and 



I = 2 equations of (56) with the Bardeen equation (81) 



we generate an equation for H u . As our gauge-invariant 
A and V are divergent, we must carefully substitute them 
out. After this we find a solution of the form 



when the growth before and after decoupling is included. 
The compensated mode is of amplitude n^'. The vector 
mode has no equivalent passive mode as perturbations 
purely to the vector potential <7 (1) decay away; it does 
have a compensated mode, again of amplitude Tig '. For 
more details see [8]. 



Statistics 



1-^ 



>s (a In (r/r v )) 
+ di sin (a In {t/t u )) 



(85) 



where a is a positive constant depending on R v . As 
t — > oo we can see that the solution Yl v —> — ^Hb, 
compensating the magnetic anisotropic stress. When the 
compensation is effective the source becomes zero and 
the potentials stop growing. The further growth in the 
curvature can be calculated giving the final curvature 



C«C(T B )-^n B 



log(Wr B ) + ( — - 1 



, (86) 



where we have neglected terms in tb/t u <C 1. 

Our initial conditions are given in the synchronous 
gauge and thus for calculations we need the curvature 
perturbation in this gauge. It can be calculated from 
Q = ?;+r)/2'H (in radiation domination) . On superhorizon 
scales, when the compensation is complete, the derivative 
term will be zero, and t](t) rs C( t )- 

At some later time when the anisotropic stress is com- 
pensated their are effectively two types of perturbation. 
The first is an adiabatic-like mode with an amplitude £ ~ 
— i? 7 II B log (t„/t b )/3, the so-called passive mode, with 
all species having zero initial anisotropic stress and un- 
perturbed densities. As we will see later, whilst the pas- 
sive mode gives adiabatic type perturbations, the statis- 
tics of TIb are non-Gaussian unlike the standard adia- 
batic mode, and will have significant higher order statis- 
tics [24] . The second type is the well known compensated 
magnetic mode (see [25 - 27 ) . with no initial curvature 
but containing the perturbed density and anisotropic 
stresses (with the total density and anisotropic stress un- 
perturbed). We consider this in two parts: an anisotropic 
stress sourced mode, with the compensating anisotropic 
stresses and unperturbed densities, and a density sourced 
mode with unperturbed anisotropic stresses but compen- 
sating densities. These have amplitudes proportional to 
lis, and Ab respectively, and their initial behaviour is 
presented in detail in Appendix [B] Though we split them 
in two, these two compensated modes are not indepen- 
dent; we address the statistics of this in the next section. 

The situation for the tensors is similar with, resulting 
in a passive tensor mode of amplitude 



H< 



log(Wr B )+ ( gj^-l 



(87) 



The statistics of Bi are assumed to be gaussian, and 
as we do not include helical fields in our analysis [25] , 
described by a power spectrum Pe(fc) defined by 



where Py = 5ij — kikj is a projection tensor that comes 
from the zero divergence of B. Calculating the energy- 
momentum perturbations requires us to consider them 
in harmonic space, and this turns the real-space mul- 
tiplications of B into fc-space convolutions. This can 
then be used to calculate the power spectra of the en- 
ergy momentum perturbations in terms of convolutions 
of the magnetic field power spectrum Pb- Various re- 
sults have been calculated for this, from approximations 
[2"T1 12U1 15U] to exact results for specific magnetic spectral 
indices pj [25]. Since the energy momentum perturba- 
tions are quadratic in the magnetic field, they cannot 
be Gaussian. Nonetheless the predicted power spectrum 
is still interesting observationally, though more informa- 
tion is available by also looking at higher-point statis- 
tics [24 l l3T | [32]. 

Though there are two scalar magnetic sources, as they 
are both sourced by the same underlying magnetic field, 
they are not independent, and when considering their 
affect on the CMB we must carefully set up the initial 
conditions for them with the correct amplitudes and cor- 
relations between them, as well as the correct relative 
amplitude of the vector and tensor contributions. 

The scalar energy density perturbation is defined 
above. The scalar anisotropic stress perturbation is 
Ii B = -|lij(fe)n^ where we denote the traceless ten- 
sor Tij(k) = (kikj — j^Sij). In terms of the magnetic field 
these are written as 

1 



A R = -&,-A 



n B = -T tj (k)A* 



(89) 



where we have hidden the convolution of the magnetic 
field in a definition of 



A lJ = 



j / d 3 pd 3 q B\p)B>(q)6(k-p- q ) 



47r(27r) 3 p 7 



< 9 °) 

There are three power spectra that we will need to com- 
pute, the power spectra of both Ab and n B , and also, 
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the oft-neglected cross correlation of the two. In terms 
of two point statistics of A 1 -? 

<A B (k)A B (k')) = ^ m (A^'(k)A im *(k')) . 

(A B (k)n B (k')> = 9 -S t] TUk'){^ (k)A^(k')) , 

(91) 

Ol 

(n B (k)n*(k')) = -T l3 {k)T lm {k') (A'^(k)A' m *(k')) . 

To calculate (A«(k)A' m *(k')) we substitute the defini- 
tion ( 90 1 , and then using Wick's theorem to evaluate the 



4-point correlator of the gaussian B, we end up with a 
result in terms of a convolution of P B 



(A ij '(k)A im *(k')) 
^( k - k ') f 

16(27r) 2 p2 a 8 

x \p a {p)P jm {q) + P m (p)P3 l (q) 



d 3 pd 3 qP B (p)P B (q)S(k - p - q) 



(92) 



With this result we can calculate the power spectra of 
the scalar perturbations by performing the relevant con- 
tractions of Tij , Pij and Sij , which leave terms dependent 

on the angles between k, q and p (or (k — q) as it will 
become when we integrate out the Dirac-delta function). 
We will denote 7 = fc ■ g, j3 = k ■ p and pi — p ■ q. The 
three correlations can be written in terms of exact inte- 
grals, first 



<A B (k)A* B (k')) = 



<J(k - k') 

1287T 2 p2 a 8 



d 3 qP B (q)P B (\k-q\)(l + n 2 ) , (93) 



second 



(n B (k)n B (k')) 



9<5(k - k') 

327T 2 /9 2 a 8 



d 3 qP B (q)P B {\k 
3 „ 1 



l-\{l 2 +p 2 ) + \l 2 p 2 -ll^+^ 2 



and lastly the cross correlation 
35(k-k') 



(94) 



A B (k)n B (k')) = 3 ^ 2 /J J d 3 q p B ( q )p B (\k-d\) 



1 



647r 2 /3 2 a : 



(95) 



In the literature, the magnetic anisotropic stress II B is of- 
ten replaced by the Lorentz force, given, in our notation, 
by L B — I (wjH B — A B /2). By combining the correla- 
tions of (A B (k)A B (k')) and (A B (k)n s (k')) we can see 
that in general there is a non-zero correlation between L B 
and A, that has often been neglected in the literature. It 
is given by 

<A B (k)L B (k')) = /^Va 8 / d3 9PB(q)PB(\k- q \) 
x [1 - 2( 7 2 + (i 2 ) + 2 7 /?/i - M 2 ] , (96) 



and should be included when calculating the effects that 
magnetic fields have on the CMB. 

We can calculate the relevant correlations for the vec- 
tor and tensor perturbations 11^ = — 6k^e[^ 1} A 1 - 7 and 

II^ = — 2e i j 2) A i: > in the same manner. The vector cor- 
relation is 



<n fl >(k)n fl >*(k')> 



18(5(k-k') 
647r 2 j o 2 a 8 



j d 3 qP B (q)P B (\k-q\ 



x [l- 2 7 2 /3 2 +7/3/i] , (97) 



and the tensor correlation is 



<nf (k)n^(k')) = 3 4^m I d 3 q P L M^A k- q > 



647r 2 p 2 a 8 

x (l+ 7 2 )(l + /3 2 ). (98) 

Our results are in agreement with those in the literature 

era Gag. 

The exact form of the magnetic power spectrum P B (k) 
is highly dependent on the production mechanism. We 
follow the rest of the literature in choosing to use a power 
law description 



P B (k) = Ak nB 



(99) 



for k < kjj, and zero otherwise. The cutoff wavenumber 
k B comes from the fact that radiation viscosity leads to 
damping of small scale magnetic fields. This is the order 
of the Silk-damping scale times the dimensionless Alfven- 
velocity [33J EH], which is time dependent, though we 
are mainly interested in perturbations sourced around 
and before recombination. The amplitude A is defined 
in terms of the expected field amplitude B\ smoothed on 
a scale A (we use the conventional A = IMpc). This gives 



.4 



(2 7 r) TlB + 5 J B 2 



r(2 



(100) 



For illustration we shall focus on nearly scale-invariant 
magnetic field spectra, since these are the only ones 
likely to give signals in the CMB on acoustic-oscillation 
scales [3D1 US] • It should be noted that it is difficult for 
causal mechanisms to give such spectra, and so to pro- 
duce large scales modes we are likely to need some infla- 
tionary mechanism. 

For scale-invariant spectra the contributions of interest 
are then from scales much larger than the damping scale 
&d, for a spectral index n B < —3/2, the cutoff becomes 
largely irrelevant. The effect on the CVs from modifying 
the power spectrum at these scales is small. For the 
compensated modes it is around 1 percent at I ~ 2000, 
and less than 3 percent at I ~ 5000. The effect on the 
passive modes will be negligible as the magnetic damping 
scale is tiny at neutrino decoupling. 

Ignoring the cutoff in the definitions of P B allows us 
factor out the /c-dependence of the above integrals and 
make them dimensionless, depending only on the spectral 
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index. For instance the integral in ( 93 1 can be rewritten 
as 



TT Mode 



J d 3 qP B (q)P B (\k- q|) (l + fi 2 ) = 2nk 



2n B +3 



du d 7 u nB (1 - 2w 7 + u 2 ) nB ' 2 (1 + ii 2 ) (101) 

10 J-l 

where we have substituted u = q/k. The angular func- 
tions [i and j3 can be written in terms of 7 and u as 

1 




7 



(1 -2u 7 + u 2 ) 1 /2 ' 

1 — "fU 

(1 -2U7 + U 2 ) 1 / 2 ' 



(102) 
(103) 



The same can be done for all the correlations above (93 1- 
(98). Whilst the integrands have singularities at u = 



(corresponding to q = 0), and u = 1, 7 = 1 (correspond- 
ing to k — q = 0), the integrals are convergent provided 
that riB > —3. We use a series expansion to integrate 
small regions around each of the poles, and numerically 
integrate the remainder. We use a nearly scale invariant 
power spectrum with ng = —2.9 giving power spectra 



PA B (k) 

PAn B (k) 
Pu B (k) 




2n B +6 



2n B +6 



9(26.30) 
3(105.55) 



2n B +6 



2n B +6 



2n B +6 



(104) 



where our power spectra are defined in a dimen- 
sionless manner (A B (k)A^(k')) = 2tt 2 (2tt) 3 S(k - 
k')k~ 3 PA B (k). The numerically calculated value is 
wrapped in parentheses. Note that these power spec- 
tra only include one of the two separate modes for the 
vector and tensor type perturbations. The shape of our 
power spectra are identical to the commonly used ap- 
proximations of [27] , but our integration predicts signifi- 
cantly different amplitudes. Using the same approxima- 
tion scheme as [37] we would predict that the angular 
integrals are equal to 2n/ (n + 3)(2n + 3). For our spec- 
tral index this is approximately ~ 20.7. Comparing this 
to the numerical results above (shown in parentheses), 
we see that the difference is up to around five times for 
the tensor power spectrum. 

In Figure[T]we show the effect that the cross-correlation 
between As and Hb has on the CMB. Despite it being 
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FIG. 1: The scalar power spectra with and without the cross- 
correlation between As and lis • Inclusion of it in calculations 
gives a consistent increase in power of around 15-25 percent 
at all scales. 



an anti-correlation we can see that it boosts power on 
all scales, as many of the perturbations are effectively 
sourced by the Lorentz force Lb = I (iv^Hb — Ag/2). 



C. Numerical Calculation 

In Figure [2] we plot the four CMB power spectra for 
primordial magnetic fields. We use the constraint of [6], 
of B\ = 4.7 nG at a scale of A = 1 Mpc, with a re- 
alistic neutrino mass m i> — 0-47 eV taken from the 
recent constraints of [Mj. We include both the compen- 
sated modes for all three perturbation types as well as 
the passive modes. Note that within this paper we as- 
sume that the magnetic perturbations are uncorrelated 
with the primary sources of anisotropy in the CMB. 

There is currently no leading theory of the formation 
of primordial magnetic fields, though there is much work 
suggesting their production could be around the elec- 
troweak phase transition [37] at T ~ 1 TeV, or from 
just after the Quark-Hadron phase transition [38J [39] at 
T ~ 150 MeV. However, to produce a scale-invariant 
spectrum we are likely to need some kind of acausal infla- 
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tionary method [10], though many often struggle to pro- 
duce large enough magnetic fields. For an unknown infla- 
tionary mechanism the exact time and details of magnetic 
field production are unclear, however, for illustration we 
believe that the electroweak transition provides a use- 
ful bound on the latest production time, and reheating 
(at temperature T < 10 14 GeV, about the GUT scale) a 
bound on the earliest. This gives t„/tb ~ 10 6 -10 12 . Any 
magnetic perturbations directly generated during infla- 
tion will source passive modes which are essentially just 
a component of the primordial spectra. 

Our CMB power spectrum results are in broad agree- 
ment with previous work |Hj in the cases where the 
results have been calculated. 

The most significant magnetic contributions to the 
CMB come from the tensor passive modes in all four 
power spectra. This is at the level of 10 percent for the 
temperature anisotropies and around an order of magni- 
tude greater than the primordial gravitational wave con- 
tribution to the B-mode polarisation. The compensated 
vector mode is important on very small scales in the tem- 
perature power spectrum [41], though here we also have 
to cope with significant secondary contributions from SZ 1 
and CMB lensing. The vector mode also leaves a clear 
signature in the B-mode polarization spectrum on small 
scales, with a comparable amplitude but different shape 
to the secondary signal expected from CMB lensing. 

Its large amplitude at low multipoles mean that the 
passive mode may provide stronger constraints on any 
primordial magnetic field than the compensated mode, 
though the relative amplitude between the two is uncer- 
tain due to the unknown epoch of magnetic field produc- 
tion (but the dependence is only logarithmic) . Using cur- 
rent WMAP temperature data, the passive modes should 
constrain the magnetic field to lower than the current 
linear-theory CMB-only limit B\ < 4.7 nG of [5]. Planck 
B-mode data will only enhance this. The effectiveness 
of these CMB constraints will be limited by: secondary 
effects at small scales obscuring the compensated vector 
mode; and confusing primordial tensor modes with the 
passive modes on large scales (with a large cosmic vari- 
ance). We should also note that as the amplitude of the 
power spectra scales like Bf, improving the upper lim- 
its on the magnetic modes, translates into much weaker 
improvements in the magnetic field strength constraints. 

The results of [9] suggested that the presence of mas- 
sive neutrinos led to a significant enhancement in power 
in the compensated modes at the largest scales. Whilst 
we see an see an increase in power on large scales due to 
massive neutrinos, the effect we calculate is much less sig- 
nificant (by about five orders of magnitude). We believe 
this effect is due to a numerical issue with tight-coupling 
which we discuss in detail in Section fVlH Al 



VIII. NUMERICAL ISSUES 

A. Tight Coupling 

To derive a tight coupling approximation for tensors, 
wc take the evolution equations for the CMB temperature 
and E-mode polarization quadrupolc 
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we rearrange these two 



equations, and substitute for E 2 to give 
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Looking at this, we see that for small kr c and small t c /t, 
the temperature quadrupole is also small, though we 
don't show it this is also true for the E-mode quadrupole. 
Physically this can be interpreted as the photons are 
tightly coupled to the baryons if there are many scatter- 
ings within a wavelength of the perturbation, and there 
are many scatterings across the horizon size. Rearrang- 
ing the equation for higher temperature moments 
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we see that higher moments are suppressed by factors of 
kr c , that is ; <2) cx kT c 9f\. If we want to only retain 
terms up to first order in r c , this allows us to ignore 



higher moments in (107). Noting that B 2 ' cx kr c Ek by 



the same argument, we can drop all terms cx kr c leaving 
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If both kr c and t c /t are small, then the -Kf an d 
terms in the right hard bracket are small corrections to 
the value of overall 0^ and we can neglect them leaving 



-T r H 
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(110) 



1 Recent work has suggested a strong constraint from the magnetic 
mode contribution to the SZ effect 1421 



the standard tight-coupling approximation for the ten- 
sors. 
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FIG. 2: The four CMB power spectra plotted for a realistic neutrino mass ~}2m v — 0.47 eV, with a magnetic field B\ — 4.7nG. 
We include the scalar primary contribution for the TT,EE and TE power spectra, and the tensor primary (with a tensor to 
scalar ratio of 0.1) and for the BB power spectrum. The shaded regions represent the regions we would expect the passive 
modes to lie within for production between the reheating and the electroweak transition. 



The problem within CAMB (Feb 2009 version) is that 
for tensor modes with small kr c it uses the tight coupling 
approximation no matter what the value of t c /t. This 
clearly invalidates the tight coupling approximation, but 
it does not manifest itself for standard models as most 
of the quantities we are interested in are proportional to 
k anyway, and thus the overall error is small (generally 
much smaller than 1 percent on the largest scale Cj's). 

As can be seen from the initial conditions in Ap- 



pendix |EJ the growth of modes like the tensor compen- 
sated magnetic mode is modified and they grow pro- 
portional to k 2 s T 2 with an effective wavenumber k 2 s — 
k 2 +arh 2 and thus on very large scales k 2 s oc m 2 . This de- 
generate evolution ensures that the growth of large scale 
perturbations remains large, and thus there is a large 
error from the tight coupling approximation. 

In Fig. [3] and Fig. [4] we show the vector and tensor 
contributions to the four CMB power spectra before and 
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FIG. 3: The compensated vector contributions to angular 
power spectra of the temperature and polarization of the 
CMB. For each spectrum we plot three different cases, for 
purely massless neutrinos (dashed), and for massive neutri- 
nos C^2m v = 1.8eV) calculated using the CAMB defaults 
(dotted), or our modified version (solid). In all cases we use 
a magnetic field of B\ = 4.7nG. We also include the primary 
contribution to the spectrum in each case (thick solid), scalar 
perturbations for the TT, EE, TE plots, and the gravitational 
wave contribution to BB. Whilst both massive neutrino cases 
contain significant large scale power compared to the mass- 
less neutrinos, our modifications avoid the artificial increase 
at very low / given by the CAMB default. 

after correcting the tight coupling behaviour. As we can 
see this has significant effects, most notably on the tensor 
contribution to the EE mode power spectrum, compared 
to the default behaviour of CAMB. This explains the 
tremendous increase in large scale i?-mode power seen 
in the results of [9] where we have used the same total 
neutrino mass = 1.8eV, magnetic field strength 

B\ = 4.7 nG and magnetic spectral index ns = —2.9. 
Our calculation shows that, even at the lowest multi- 
poles, the tensor compensated magnetic mode is signifi- 
cantly lower amplitude than the primary scalar adiabatic 
spectrum, and remain subdominant to the compensated 
vector mode. 



B. Early-Time Numerical Instabilities 

In Fig. [5] we show the evolution of the tensor pertur- 
bations of several different large scales for the compen- 
sated magnetic mode. The top set of panels show the be- 
haviour in the presence of massless neutrinos, and we see 



TT Mode 4 EE Mode 




FIG. 4: The compensated tensor mode to the four CMB an- 
gular power spectra of temperature and polarization. This is 
the tensor equivalent of Fig. [3] The solid line is our modified 
version, dotted the CAMB default and dashed the massless 
case. The CAMB default exhibits the same small / excess as 
in the vector case, and as before our modified version avoids 
this. 



the slowly growing scale-dependent evolution of both the 
gravitational waves and the total anisotropic stress. The 
middle panels show the output of CAMB when evolv- 
ing massive neutrinos, illustrating a fundamental prob- 
lem when numerically evolving these perturbations. To 
evolve neutrino quantities such as the anisotropic stress, 
we need to evolve the distribution function perturbation 
^v(ki,q,nj) at a fixed set of points q, then we numer- 
ically integrate over the points to calculate the desired 
quantity. For the standard modes this approach is fine, 
however, in the case of the compensated magnetic mode, 
the initial cancellation is at the order of 10 -10 , and re- 
quire numerical accuracy at this level to calculate the 
anisotropic stress correctly. As well as simple numerical 
accuracy, we must integrate well into the tail of the dis- 
tribution to include all contributions up to 10 -10 , this 
requires an increase in the range of q values integrated 
from g max ~ lbkgT up to around g max ~ 40fcsT. 

One way to obtain the numerical accuracy would be to 
simply increase the number of points over which we inte- 
grate. However, combined with the required increase in 
range, this requires a significant increase in the number of 
integration points. We use an alternative approach, using 
our mass expansion of the Boltzmann hierarchy. At early 
times we directly evolve the integrated moments Fi and 
use this to calculate the anisotropic stress. As the neutri- 
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FIG. 5: The evolution of the tensor metric perturbation i/ <2) 
(left panels), and the total anisotropic stress II <2) plotted 
against the scale factor a, at various wavenumbers. In the top 
panel we show the evolution with massless neutrinos. The 
middle panels illustrate the behaviour when we instead use 
three massive neutrinos X^ m " = 0.18eV, with the default 
behaviour of CAMB. The problems stemming from the in- 
tegration accuracy are readily apparent at early times. The 
bottom panels show the correct evolution of the massive neu- 
trinos with our modifications. The degenerate evolution at 
small k is apparent. 



nos start to become non-relativistic, our mass expansion 
becomes inaccurate and so before this we switch to using 
the full distribution function. By this time the level of 
cancellation is within the numerical accuracy of the inte- 
gration and the total anisotropic stress is accurate. The 
results of this are shown in Fig. [5] 

Such an approach is essential to accurately model the 
behaviour of massive neutrinos in the early Universe, 
however for calculating CMB power spectra the correc- 
tons are sub-percent level and simply increasing the range 
and number of integration points is sufficient. 



IX. CONCLUSION 



In this paper we have developed an integrated Boltz- 
mann hierarchy for analysing massive neutrinos in the 
early Universe which is accurate to second order in the 
mass. We have calculated the leading order mass cor- 
rections to the initial series solutions for the regular per- 
turbation modes, and also demonstrated its use for accu- 
rately evolving massive neutrinos in the early Universe. 

We have made a detailed analysis of the effects of the 
primordial magnetic fields on the CMB. In our examina- 
tion of the statistics of the magnetic field perturbations 
we have included an often neglected cross-correlation 
term between the two scalar-perturbations. This serves 
to increase the power in the CMB from the compensated 
mode by around 25 percent at all scales. We also demon- 
strate that one of the standard approximations to the 
statistics can give an amplitude around a factor of five 
smaller than a more accurate result, reinforcing the need 
to move to more advanced results. 

We accurately calculate the contributions of the vari- 
ous magnetic modes (both passive and compensated to 
the CMB). By correcting some numerical issues we come 
to different conclusions to jjj] . Whilst we agree that there 
is an enhancement to the large scale power spectra (espe- 
cially E-mode polarization) caused by massive neutrinos, 
we find a much smaller amplitude; too small to enhance 
prospects of detecting primordial magnetic fields. Our 
work suggests that the passive modes are likely to dom- 
inate the compensated modes with a power spectrum 
amplitude several orders of magnitude greater at large 
scales. With the magnetic field we have used they are 
around 10 percent of the primary spectra, and this sug- 
gests that they will provide the biggest constraint on any 
primordial magnetic field in the near future, adding in 
a small gravitational wave like component with a blue 
spectral index. However unlike the compensated mode, 
such modes are dependent on the details of the mag- 
netic field production, though quite weakly, and cannot 
provide model independent constraints on the magnetic 
field in the same manner as the compensated modes. 

The Mathematica notebooks to calculate the mag- 
netic mode amplitudes, and the initial conditions can be 
found at http://camb.info/jrs/. Our modifications to 
CAMB to calculate the compensated magnetic modes are 
at the same location. 
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Appendix A: Higher-order mass expansion 



Wc define a scaled mass 



(Al) 



so the ratio of massive and massless neutrino densities is 
given by 

Pu 120 



where 
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Performing an expansion of /(to) in the mass by per- 
forming a series expansion of the square root inside the 
integral is not valid since to is not much smaller than q 
over the full range of the integral. Instead we split up 
the integral at a point a (where to <C a -C 1) so that 
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The result is independent of a, and evaluates numerically 
to 
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Thus the next term above the leading mass correction we 
consider in the paper is 0(fh A ln(rn)). A similar approach 
can be followed for a mass expansion of the pressure using 
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Appendix B: Initial Conditions 

Here we present initial series solutions for the regu- 
lar modes for scalar, vector and tensor perturbations in 
cosmology. We allow for two significantly different neu- 
trino mass eigenstates, allowing us to describe most of 
the possibilities of the neutrino mass hierarchy. The so- 
lutions are correct to order to 2 in the neutrino mass. For 
space we have only included the terms up to second order 
or the first non-zero term up to order r 3 . 

We include the standard matter species which we gen- 
erally denote with subscripts: photons (7), baryons (b), 
cold dark matter (c), massless neutrinos (v) and massive 
neutrinos (n). Our solutions are for after neutrino de- 
coupling; we discuss the pre-decoupling behaviour in the 
presence of magnetic fields in the Section [Vll A| To solve 
the evolution of the background equation we solve the 
Fricdmann equations for the scale factor. The solution 
to order to 2 in the neutrino mass is 



<z(t) = a 



1 1 2 2 
LOT + -10 T 



1 n 2 

— u r 
96 "ft 2 



- 2 4 4 

TO LO T 



(Bl) 



where we choose some time to fix the values of tt r — 
Q~ + Slj, + ft rl and f2 m = flj, + Q c and R n — Q n /il r . 
We have used the standard definition of Q x = p x /p C i the 
ratio of the density of species x to the critical density. In 
the above we also use the definition of to — r2 m 'Ho/v / ftr- 

To give our series solutions we will also use the defini- 
tions of R 1 = fi 7 /f2 r , R v = Qj,/n r , R t — (fl v + f2„)/Q r , 
R c = fi c /Q m and R b = Vl b /Vt m . 

The solutions were calculated using Mathematica, and 
the notebooks and relevent packages used can be found 
at |http: //camb. inf o/jrs| 
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1. Scalar Initial Conditions 



There are six regular scalar modes, one adiabatic, four isocurvature, and one magnetic. For comparison to other 
results we give our solutions in the synchronous gauge [3143] with the standard potentials h and rj, commonly used for 
its numerical robustness. This has the further advantage that the neutrino velocity isocurvature mode is completely 
regular as r — > [33]. We also give the Bardeen potentials used in the text and 
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Baryon Isocurvature Mode 



Baryon isocurvature modes are essentially observationally indistinguishable from a rescaled CDM isocurvature 
mode. This is because the compensated mode (with 6p b = —Sp c ) gives only a small contributions at small scales, 
primarily from the baryon pressure and second order effects |44) . 
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Neutrino Isocurvature Mode 



3R ^ 

^ t) = {-w^r7) — mr) T +0{T) 

W=i + (4- (2 + ^ 2 ^ )^ + o ( ,3) 

-n(r)-^+0(r 3 ) 

QjL2 2 
1.3 3 

MT) ^ i+(4 _^p )r2 + 0(T3) 

*v(t)^+0(t 3 ) 

1.3 3 

/fc 2 ^ 3i?„c 2 m 2 2 ^ 2 3 
, , (fci? t )r , ikR b R t ujT 2 3 

" b(T) ^'^7 + i6A 2 +0(T) 



21 



2 



%(t) _<^ + 3^ + 0(tS) 

t/ n 2i? t (75-2R t )R t LJT _ ,, 

*fr) = 1 — h O t ) 

[ ' 15 + AR t A (225 + 90i? t + 8i? t 2 ) ^ v ; 

*(r) = flt + ^(-15 + 2^ 2 
W 15 + ARt A (225 + 90R t + 8R 2 t ) { ' 



Neutrino Velocity Isocurvature Mode 

Despite the apparent singularities in the potentials <J/ and <&, the mode is physical with a regular comoving curvature 
perturbation. However, as the neutrinos are strongly coupled to the photons prior to decoupling it is challenging to 
find a mechanism to source this mode. 
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4flt flf(-15 + 4fit)h; 

lTj fc(5 + 4iZ t )r + fc(5 + 4J2t)(15 + 4i? t ) + U(Tj 

In theory we can define isocurvature modes in the neutrino anisotropic stress and higher multipole moments, though 
with no reasonable mechanism to produce them we will omit them. 



Compensated Magnetic Modes 



For the compensated magnetic mode, we treat it like an isocurvature mode with 77 — > at very early times. For 
the density sourced modes this gives 
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The anisotropic stress sourced modes are 
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2. Vector Initial Conditions 

There are two regular vector modes, a vorticity mode which is the vector equivalent of the neutrino velocity 
isocurvature mode, and a magnetic mode compensating the magnetic anisotropic stress H ( g . We give the solutions 
in terms of the gauge invariant variables used earlier. 



Vorticity Mode 

For the same reasons as the neutrino velocity mode, the existence of this type of perturbation is highly unlikely. 
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3. Tensor Initial Conditions 

Only the photons and neutrinos can support tensor perturbations to their energy momentum tensors and at times 
long before recombination the photon anisotropic stress is negligible. Thus the species affecting the tensor evolution 
are the neutrinos, and the magnetic fields. This leaves us with one standard tensor mode, the gravitational wave 
mode, and a compensated magnetic mode. 



Gravitational Wave Mode 
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